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We have created a fiat piling of disks in a numerical experiment using the Distinct 
Element Method (DEM) by depositing them under gravity. In the resulting pile, we 
then measured increments in stress and strain that were associated with a small decrease 
in gravity. We first describe the stress in terms of the strain using isotropic elasticity 
theory. Then, from a micro-mechanical view point, we calculate the relation between 
the stress and strain using the mean strain assumption. We compare the predicted 
values of Young's modulus and Poisson's ratio with those that were measured in the 
numerical experiment. 
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1. Introduction 

The ultimate goal of our research is to describe how stresses propagate through granular 
media. Physical experiments 1 ) indicate that when a localized force is applied at a surface of a 
layer of sand, the stress distribution on the bottom of the layer depends on how the layer was 
constructed. The pressure distribution under a conical sand pile also depends on the way the 
pile was built and may have a local minimum under the apex of a pile. 2,3 ) We don't yet know 
how to describe the way that stress propagates through a pile or how the propagation depends 
on its construction history. 

In this paper, as a first step, we focus on a two-dimensional example with possible anisotropy 
in the vertical direction: a flat pile of circular, elastic, frictional disks, deposited under gravity 
onto a frictional, flat base. In this case, as shown in § 4.3, we find that the anisotropy due to the 
gravity is small enough to neglect. In § 2, we write down equations of force balance, assuming 
that the pile is continuous body, in order to predict stress and strain within it. Assuming that 
the pile is an isotropic elastic body, these equations are solved exactly. But the problem is that 
material constants, such as the Young's modulus and Poisson's ratio are unknown. In the next 
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chapter, we adopt the mean strain assumption in order to predict the material constants of 
an assembly of disks in terms of micro-scopic material properties, such as the contact stiffness. 
Here, a relationship between the micro-scopic and macro-scopic material properties is obtained. 
In § 4 we build a flat pile with elastic, frictional circular disks in a numerical experiment using 
the Distinct Element Method, a kind of molecular dynamics. Then, the stress and strain in 
the pile are measured in a deformation produced by reducing gravity, and the predictions that 
are made in § 3 are tested against the measured values. In the last chapter, we discuss the 
reasons for the differences between the predictions based on the mean strain assumption and 
the measured results and indicate what more should be done. 



2. Elasticity theory 

2.1 Continuum theory 

The equation of motion for a granular piling is written as 

do~ij 

-^f- + pbi = p cii, 

where aij is the stress, p is the mass density, bi is the body force per unit area, and is the 
acceleration. In this paper, we focus on only a two-dimensional case, so the summation should 
be taken from 1 to 2. The body force is gravity having only a vertical component, b 2 = — 9- 

We are interested in a static piling, so a« = 0, i = 1,2. The equilibrium equations for the 
two-dimensional case are 

pl + d -p^ = 0, (1) 

OX\ 0x2 

da 12 da 22 , n . 
-dxT + ^- p9 = °" (2) 

The boundary conditions are defined for a pile of uniform height h: the horizontal boundary 

is periodic and, at the free surface, the shear stresses and the vertical stress are zero, 022(2:2 = 

h) = and 012(^2 = h) = 0. Then, by symmetry a\ 2 = 0. By solving eqs. (1) and (2), the 

components of stresses are 

O-ll = (Tll(x2), (3) 

022 = P g{x 2 - h). (4) 

In these equations, 011 is indeterminate, so the exact solutions can not be obtained using only 
this information. In the next section, we will adopt isotropic elasticity theory in order to obtain 
an exact determination of a\\. 
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2.2 Isotropy 

The horizontal and vertical strains, En and E22, can be obtained by introducing the con- 
stitutive relations of isotropic elasticity: 

En = -^(oii - 1/0-22), ( 5 ) 
E22 = -^{(?22 ~ van), (6) 

where v is Poisson's ratio and E is Young's modulus. 

The horizontal stress, an, given by eq. (3), is obtained using eq. (5) and requiring that 
En = because of the horizontal periodic boundary condition. Then 

cr n = VO22 

= vp g(x 2 - h). (7) 

Substituting eqs (4) and (7) into eq. (6), the vertical strain is 

1-v 2 

E>22 — ^ C22 

= ] —^ P g{x2-h). (8) 

3. Prediction of the material constants 

From a micro- mechanical view point, the average stress can be written using an orientational 
distribution of contacts, D(n), and a contact force, f c , as 



an = — / D(n)f?n j( m, 

7TO J 



where 7 is the area fraction of the disks, a is a radius of disk, n is the unit vector from the center 
of the disk to a contact, and dQ is the element of contact angle. 

Here, we make the strong assumption that the contact displacement is determined by the 
average strain as 

Up — (lEpgTlq. 

Then, the contact force is also determined by the average strain as 

fi = KipUp 

— E^ipdEpqTlq, 
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where 

K ip = K N mn p + K T (Si p - riiUp), 

in which Kjy and Kt are normal and tangential stiffness, respectively. 

We assume that the orientational distribution of contact is isotropic; then 

D(n) P) 

where A; is the average number of contacts per disk. 10 ) As shown in § 4.3, it is appropriate to 
assume that the contact angle distribution is isotropic. 
Then, 

cry = — / —aK ip n q njdQE pq . (10) 
7ra y Z7r 

The general form of Hooke's law is 
By comparing eqs (10) and (11), 

7 f k 
ira J 2ir 

After the integration in eq. (11) is carried out, eq. (10) becomes 

When the stress is expressed in terms of strain, the coefficients are called Lame's constants. 
Now, Lame's constants are obtained in terms of stiffness of disks as 

2fi = 2^(K N + K T ), 

A = l!L{K N -K T ). 

Then Young's modulus, E and Poisson's ratio, is, are converted from Lame's constants, so they 
can also be expressed in terms of stiffness of disks as 

A K N -K T 



v = 



2/i + A 3K N + K T ' 



+ A) _ jk (K N + K T )K N 
2/x + A vr 3K N + K T ' 1 ' 

Using this relation between micro-scopic and macro-scopic material property, we can predict 
Young's modulus and Poisson's ratio. 
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4. Numerical experiments 



In this paper, we adopt Distinct Element Method (DEM), which was invented by Cundall, 4 ) 



in order to produce deformed granular aggregates in which stress and strain can be measured. 
We follow the algorithm of A. Shimosaka 5 ** and H. Hayakawa. 6 ** 

Using DEM, we make a two-dimensional flat pile that is composed of elastic, frictional, 
circular disks on a frictional smooth bottom, and measure stress and strain at each point. 

4-1 Setting 

Elastic, frictional, circular disks are deposited on a flat frictional bottom, layer by layer. 
Two diameters of disks are chosen that are slightly different from each other, in order to avoid a 
crystal structure. The horizontal boundary is periodic. We let this pile relax until the disks come 
to a static state, after dissipating all of their kinetic energy in collisions. (See the left panel of 
Fig. 1.) There is the possibility of the existence of depositional anisotropy in the vertical direction 
that will be considered later in § 4.3. The right panel of Fig. 1 shows a Voronoi tessellation of 
this pile. Voronoi tessellation divides the whole domain into cells using perpendicular bisectors 
of the lines of centers, so that each cell contains one center, indicated by a dot. Voronoi cells 
are used to measure strain in a granular assemblies as defined in § 4.2. 

The parameters that are used in the actual simulation are normalized so that the maximum 
disk diameter, the gravitational acceleration, and the mass per unit area are all unity. Conse- 
quently, upon taking the summation over all contacts B with disk A, the dimensionless equation 
of motion for disk A is 



which are linearly proportional to the relative displacement and the relative velocity, respectively. 
A dash denotes a non-dimensional variable. For example, a dimensionless elastic coefficient and 
a dimensionless viscous coefficient are calculated as k' = {l/mg)k and rj' = (l/fh)y/l/g 77, 
respectively, where in is the mass per unit area and m' A is the dimensionless mass of disk A. 
The dimensional parameters used in our calculations are shown in Table I. 

4-2 Measurement of stress and strain 

There exist several definitions for stress and strain of granular aggregates. 7 ~ 9 ) In this pa- 
per we adopt a simplified form of definitions, 7 ' 8 ) that are consistent with each other in two 





contact forces consist of elastic and viscous forces 
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dimensions. After reducing gravity by 10%, we measure the displacements of the disks and the 
increments of the contact forces and, from them, calculate increments of stress and strain. The 
increment of the strain at a point at the center of disk A is taken to be 

where 7 is the area fraction of the diska, u AB is a displacement of disk A relative to disk B, 
n AB is the unit vector in the direction of contact between disk A and B, and l AB is the length 
of the side of a Voronoi cell which is shared by cell A and cell B, with the name of a cell the 
same as the name of a disk which is inside of the cell. The corresponding increment in stress is 
taken to be 

• A _ _1_ ST f AC „AC 

where f c AC is a contact force exerted on disk A by disk C and cia is the radius of disk A. In 
the definition of strain, the summation is taken over all the neighbourings which share sides of 
Voronoi cells in between; while for stress, the sum is taken over all of the pairs of contacts. 

The increments of stress and strain are measured at each disk center according to the 
definitions above. They are taken to be an average over disks which are included in a horizontal 
slice with a width of 1.5 dimensionless units at each height. Figure 2 shows the increment of 
the dimensionless stress versus, the height, normalized by the maximum diameter. It is linearly 
proportional to the height of the pile due to the weight of the disks. 

Now we are interested in the increments of stress and strain associated with reduction 
of gravity by 10%. Then the increments of stress components are expressed in terms of the 
increment of gravity, 5g = — 0.1<? as 

oil = v a'22 (13) 
= vp Sg(x 2 ~ h) (14) 

and 

0-22 = p Sg(x 2 - h). (15) 

The slope of the line along a'22 coincides well with the value, —0.07, which is estimated from 
eq.(15), where p = 0.83 which is calculated from the area fraction. In this respect, the continuous 
description pictures the stress well. 
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Figure 3 shows the increment of strain versus, height. From eq.(8), the increment of strain 

is 

\-v 2 . 

E22 = — ^ — <y~22 (16) 
hi 

= l -^p5g{x 2 -h). (17) 

Comparing eq. (17) and Fig. 3, we found that £"22 is related linearly with the height of the pile 
with some fluctuation. 

Next, let us consider the relation between stress and strain. Figure 4 shows the increment 
of the dimensionless stress versus the increment of strain. From eq. (16), stress is expected to 
be linearly proportional to strain. The linear relation can be seen in Fig. 4 between stress and 
strain in our numerical experiment. 



4-3 Contact angle distribution 

Figure 5 shows the contact angle distribution. There are peaks around 7r/3, 2ir/3, and n. It 
indicates that the disk configuration is nearly a hexagonal packing. Although we expected that 
there would exist a depositional anisotropy due to the process of construction of the pile, it was 
found that the contact angle distribution can be assumed to be isotropic in a macroscopic sense. 
So the assumption that was made in eq.(9) is supported. The hexagonal structure also can be 
seen in the Voronoi tessellation (see the right panel of the Fig. 1, although it should be noted 
that a side of Voronoi cells is defined as a bisection of the nodes connected the centers of a pair 
of disks, but it doesn't mean that they are necessarily in contact as shown by the fact that the 
average number of contacts is about 4.7. 



4-4 Measurement of the material constants 

Under the assumption that the material is isotropic, we can characterize its elastic response 
using only two constants (e.g., the Lame constants, or Young's modulus and Poisson's ratio). 
However, we don't yet know such constants for the granular aggregate as a bulk. In other 
situations involving anisotropy, two constants are not be sufficient to characterize the elasticity 
of a granular material. 

For the isotropic material, we can estimate Young's modulus and Poisson's ratio from eq. (13) 
and (17) as 

o-'n 

v = —, 

0~22 
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E = (18) 

The ratio of an against 022 is almost constant as shown in Fig. 6 except at the surface of the 
pile. Also, Young's modulus, estimated from eq. (18), is shown in Fig. 5. The fluctuation in 
Fig. 5 is found is larger than that for Poisson's ratio and may not be negligible. 

As predicted in § 3, using 7 = 0.83 and k = 4.7, as measured in the numerical simulation, 
the dimensionless Young's modulus E' = (l/fhg)E is 9.3 x 10 3 and the Poisson's ratio is 0.25. 
We can compare those values with those obtained in the numerical experiments shown in Figs. 6 
and 5. The predicted Young's modulus is almost 2.58 times of the measured value; while the 
predicted Poisson's ratio is half of the measured value. That is, the predicted properties of the 
pile are stiffer than those measured. 

4-5 Fluctuation of strain 

In our micro- mechanical view point, fluctuations in the strain were not taken into account. 
Consequently, the predicted material constants were larger than those measured in the numerical 
simulation, because the disks in the simulation can translate and rotate in ways different from 
that predicted by the average strain in order to satisfy force and moment equilibrium. These 
additional degrees of freedom allows them to behave more flexibly than predicted. Figure 8 shows 
the fluctuation of strain (AE22/E22) versus the height. The fluctuation, AE22, is evaluated by 
the ratio of the absolute value of the deviation from the mean strain to the local value of the 
mean strain. On average, the strain fluctuation 5E22 is 65% of the strain itself, except near 
the bottom and the surface of the pile. The fluctuation of strain is very large, especially near 
the boundaries. We need to take this fluctuation into account in order to better describe the 
relation between stress and strain. 11 ) 

5. Conclusion 

In this paper, we have investigated a flat piling of disks in order to define the relation be- 
tween stress and strain using a numerical experiment and a continuum theory. As a continuum 
description, we adopted isotropic elasticity theory. In a micro-mechanical approach, we calcu- 
lated the stress using the mean strain assumption and predicted the Young's modulus and the 
Poisson's ratio of the pile. The material constants which are measured in the numerical exper- 
iments are smaller than those which were predicted using mean strain assumption. In order to 
obtain a better description, we need to introduce the fluctuation of stress and strain in the pile. 
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The disorderness of the disk configuration may cause this reduction of stiffness as an aggregate 
as we can make sure in a pile of disks with a square lattice structure. As a next step, we can 
introduce an anisotropy of a pile, and see how stress propagation depends on it. 

An important questions that should be addressed is the extent to which elasticity theory 
can be applied to describe unloading of the pile. Preliminary numerical experiments indicate 
that particle sliding gives rise to irreversible behavior for decreases in gravity larger than ten 
per cent. The characterization and description of this inelastic behavior will be the subject of a 
subsequent paper. 
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Table I. Parameters in our numerical experiments 



Diameters of disks 


7.6 and 8.0 (10~ 3 m) 


Thickness of disk 


6.0 (10~ 3 m) 


Density of disk 


1.06 x 10 3 (kg/m 3 ) 


Gravitational Accel. 


10 (m/s 2 ) 


Restitution Coef. 


0.6 


Frictional Coef. 


0.4 


Normal elastic Coef.(fc n ) 


1.27 x 10 4 (N/m) 


Tang, elastic Coef. (kt) 


2.54 x 10 3 (N/m) 


Normal viscous Coef.(ry n ) 


1.00 (kg/s) 


Tang, viscous Coef.(?7t) 


1.00 (kg/s) 


Time step (dt) 


3.16 x 10~ 7 (s) 
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Fig. l. 

Disk Configuration Voronoi Tesselation 
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Fig. 2. 
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Fig. 3. 

INCREMENTAL STRAIN 
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Fig. 4. 
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Fig. 6. 

Poissons ratio 
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Fig. 7. 
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Fig. 8. 

Fluctuation of Strain 
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Figure 1. Left panel: The configuration of disks consisting from 495 disks in a static state. 
The horizontal boundary is periodic. The diameters of disk is 0.76 cm and 0.8 cm. Right panel: 
The lines show Voronoi cells and dots show centers of disks. 

Figure 2. The increment of the dimensionless stress versus the height normalized by the 
maximum diameter (D), The inclinations of the straight lines for a' u and a 22 are —0.0664 and 
—0.033 respectively. This dimensionless stress, a' , is converted to the dimensional value, a, as 
a = {fhg/l) a'. 

Figure 3. The increment of strain versus height normalized by the maximum diameter (D). 
The slopes of the straight lines for En and E22 are and —1.31 * 10 -5 respectively. 

Figure 4. The increments of stress, a xx and a yy , versus the increment of strain, E yy . The 
relation between the increments of stress and strain seems to be linear. The slope of the straight 
line is 2.2 * 10 3 (for a' u (x),and 4.85 * 10 3 for a' 22 (+)). 

Figure 5. Probability distribution function of contacts versus contact angle from to ir. 
Three peaks are found 7r/3,27r/3 and n. It implies that the packing can be assumed to be almost 
hexagonal. 

Figure 6. Poisson's ratio, versus height normalized by the maximum diameter. This is the 
ratio of an against 022 The mean value is about 0.5. 

Figure 7. The dimensionless Young's modulus, E' , versus height normalized by the max- 
imum diameter.lt is estimated from eq. (18). The mean value is around 3.6 * 10 3 , and the 
dimensional value, E, is converted with the relation E = (fhg/l) E' . 



20 Shio Inagaki and James T. Jenkins 

Figure 8. The ratio of the fluctuation of E22, AE22, to E22 versus height normalized by the 
maximum diameter. The strain fluctuation is almost 74 % of the strain in average except near 
the surface and the bottom. 



